NEUTRON STAR POPULATION DYNAMICS.II: 
3D SPACE VELOCITIES OF YOUNG PULSARS 



J. M. Cordes 

Astronomy Department and NAIC, Cornell University 
David F. Chernoff 

C — ■ 

q\ ■ Astronomy Department, Cornell University 

G\ ' 

ABSTRACT 

■ We use astrometric, distance and spindown data on pulsars to: (1) estimate 

three-dimensional velocity components, birth distances from the galactic plane, 
and ages of individual objects; (2) determine the distribution of space velocities 
and the scale height of pulsar progenitors; (3) test spindown laws for pulsars; (4) 
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test for correlations between space velocities and other pulsar parameters; and 



(5) place empirical requirements on mechanisms than can produce high velocity 



o 
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neutron stars. Our approach incorporates measurement errors, uncertainties 

Q\ ■ 

■ in distances, deceleration in the Galactic potential, and differential galactic 

rotation. We focus on a sample of proper motion measurements of young pulsars 
(< 10 Myr) whose trajectories may be accurately and simply modeled. This 
+l5 . sample of 47 pulsars excludes millisecond pulsars and other objects that may 

have undergone accretion-driven spinup. 

We estimate velocity components and birth z distance on a case-by-case basis 
assuming the actual age equals the conventional spindown age for a breaking 
index n=3, no torque decay, and birth periods much shorter than present day 
periods. Every sample member could have originated within 0.3 kpc of the 
galactic plane while still having reasonable present day peculiar radial velocities. 

For the 47-object sample, the scale height of the progenitors is ~ 0.13 
kpc and the three dimensional velocities are distributed in two components 
with characteristic speeds of 1 75^30 km s _1 and 700ti5o km s _1 , representing 
~ 83% and ~ 17% of the population, respectively. The sample velocities are 
inconsistent with a single component Gaussian model, well-described by a two 
component Gaussian model but do not require models of additional complexity. 
From the best-fit distribution, we estimate that about 20% of the known pulsars 
will escape the Galaxy, assuming an escape speed of 500 km s _1 . The best-fit, 
dual-component model, if augmented by an additional, low velocity (< 50 km 
s _1 ) component tolerates, at most, only a small extra contribution in number, 
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less than 5%. The best three component models do not show a preference for 
filling in the probability distribution at speeds intermediate to 175 and 700 
km s _1 but are nearly degenerate with the best two component models. We 
estimate that the high velocity tail (> 1000 km s _1 ) may be underrepresented 
by a factor ~ 2.3 owing to selection effects in pulsar surveys. 

The estimates of scale height and velocity parameters are insensitive to the 
explicit relation of chronological and spindown ages. A further analysis starting 
from our inferred velocity distribution allows us to test spindown laws and 
age estimates. There exist comparably good descriptions of the data involving 
different combinations of braking index and torque decay timescale. We find a 
braking index of 2.5 is favored if torque decay occurs on a time scale of ~ 3 Myr, 
while braking indices between 3.5 and 6 are preferred if there is no torque decay. 
For the sample as a whole, the most-probable chronological ages are typically 
smaller than conventional spindown ages by factors as large as two. 

We have also searched for correlations between three-dimensional speeds 
of individual pulsars and combinations of spin period and period derivative. 
None appears to be significant. We argue that correlations identified previously 
between velocity and (apparent) magnetic moment reflect the different 
evolutionary paths taken by young, isolated (nonbinary), high- field pulsars 
and older, low-field pulsars that have undergone accretion-driven spinup. We 
conclude that any such correlation measures differences in spin and velocity 
selection in the evolution of the two populations and is not a measure of 
processes taking place in the core collapse that produces neutron stars in the 
first place. 

We assess mechanisms for producing high-velocity neutron stars, including 
disruption of binary systems by symmetric supernovae and neutrino, baryonic or 
electromagnetic rocket effects during or shortly after the supernova. The largest 
velocities seen (~ 1600 km s _1 ) along with the paucity of low-velocity pulsars 
suggest that disruption of binaries by symmetric explosions is insufficient. 
Rocket effects appear to be a necessary and general phenomenon. The required 
kick amplitudes and the absence of a magnetic-field-velocity correlation do 
not yet rule out any of the rocket models. However, the required amplitudes 
suggest that the core collapse process in a supernova is highly dynamic and 
aspherical and that the impulse delivered to the neutron star is larger than 
existing simulations of core collapse have achieved. 

Subject headings: pulsars, stars-binary: 
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1. INTRODUCTION 

The high velocity nature of pulsars has been recognized since the first, minimal 
sample was analyzed by Gunn & Ostriker (1970). As pulsar discoveries and proper-motion 
measurements have increased, and the pulsar distance scale improved, the largest inferred 
velocity has increased to ^1000 km s _1 . The population's dispersion exceeds that of any 
other normal Galactic stellar population. 

Although it is likely that a pulsar's characteristic motion is generated during the 
neutron star formation, a number of physical mechanisms has been proposed. In broad 
terms, these include: (1) disruption of binaries through instantaneous, symmetric mass 
loss in supernova explosions (Blaauw 1961; Gott, Gunn & Ostriker 1970; Radhakrishnan 
Sz Shukre 1985); (2) slow, post-explosion, rocket effects associated with loss of spin energy 
by the neutron stars (e.g. Harrison & Tademaru 1975; Helfand & Tademaru 1977); and 
(3) an instantaneous momentum impulse, or "kick", imparted through asymmetry of the 
supernova explosion (Shklovskii 1970; Dewey & Cordes 1987). 

Many statistical analyses of pulsar velocities have attempted to discern which (if 
any) of the mechanisms is favored (Dewey & Cordes 1987; Bailes 1989; Iben & Tutukov 
1996). The recent treatment by Lyne & Lorimer (1994; hereafter LL) has incorporated new 
proper-motion measurements (Harrison et al. 1993) and an up-to-date model for the free 
electron density in the Galaxy (Taylor & Cordes 1993; hereafter TC). The latter is essential 
for estimating distances to most pulsars. LL find a larger mean pulsar speed, V ~ 450 km 
s^ 1 , than had been estimated previously. The probability density function (pdf) for the 
perpendicular speed is 

04AV - 3 
i + (vjv 

where Vb = 330 km s~ x (and dN oc fv ± d 2 V±) []. LL state that ~ 50% of pulsars will escape 
the Galaxy, assuming an escape speed ~ 500 km s _1 . LL also re-emphasized the fact that 
pulsar surveys tend to undercount the highest velocity pulsars. Surveys are signal-to-noise 
limited and sample a finite volume which, of course, is roughly centered on the Galactic 
plane. Since it is generally thought that most pulsars are born in the vicinity of the Galactic 
plane and move away from it, fast objects spend less time in the detectable volume than 
slow ones and are consequently under-represented. 

In attempting to characterize the pulsar velocity distribution, a fundamental difficulty 
is that all of the basic components necessary for a kinematic description of the population 
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1 The given expression corrects a typographical error in LL, according to Mollerbach & Roulet (1997) 
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are poorly known. Except for a few historical supernovae, the chronological age of each 
individual object is unknown and it must be estimated from the "spindown" age. The 
distance is uncertain, though it may be constrained by a pulsar's dispersion measure (DM) 
from knowledge of the electron density within the Galaxy. The perpendicular components 
of the velocity may be derived from the proper motion and distance, though the proper 
motion measurements only exist for a subset of the entire pulsar population and both 
distance and proper motion determinations have significant errors. In summary, there are 
significant uncertainties in the age, distance and velocity of nearly every known pulsar. 

Our general goal is to infer properties, like the distribution of kick velocities at birth, 
for the pulsar population. Since there is no a priori way to separate the uncertainties we 
must strive to treat them in an even-handed, consistent manner. Although we do not 
completely fulfill that goal in this paper, we make significant progress by incorporating the 
uncertainties in position and velocity directly into our analysis and by separately considering 
the sensitivity of the results to the relation between spindown and chronological age. Our 
major omission, which we plan to treat in the future, is the observational selection effects 
associated with the detection problem itself. In effect, the observer's selection function 
involves kinematic, geometric and detection-related pieces and we omit the last, focusing 
purely on the kinematic and geometric aspects. 

In this paper we present a new analysis of pulsar velocities that (1) tests and exploits the 
hypothesis that pulsars are born in or near the plane of the Galaxy; (2) takes into account 
the uncertainties in pulsar distances and proper-motion measurements; (3) estimates the 
full three dimensional velocity pdf and birth scale height of the known pulsars; (4) uses the 
global information to sharpen the estimate of measured kinematic quantities and to infer 
the most likely values of other unmeasured quantities for each pulsar; and (5) constrains 
the braking index and torque decay time for pulsar spindown. Our methodology uses a 
likelihood analysis, augmented by a comparison of models with different complexity using 
Bayesean odds ratios. This paper is the second paper in a series, where the first addressed 
the population of millisecond pulsars, finding them to be a low-velocity population with 
rms velocity ~ 84 km s" 1 (Cordes & Chernoff 1997; hereafter Paper I). 

The first half of the paper deals directly with the statistical procedures we have 
developed and the application to the data. In we describe the selection criteria for 
the pulsar sample. ^ outlines and applies our procedure for constraining radial velocities 
and birth altitudes for individual objects. §| derives a best fit three dimensional velocity 
distribution for the population and discuss the kinematic constraints on the braking index. 
§|)-0 use the statistical results as an aid in estimating separate velocity components, birth 
z and ages for individual pulsars. Using the estimated velocities we discuss in §|| the 
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velocity-magnetic moment correlation. §^ relates recent work on pulsar-supernova remnant 
associations to our analysis of the proper motion sample. 

The second half of the paper relates our work to a number of topics of current 



astrophysical interest. In § |10| we review the lines of evidence for asymmetric supernova 
kicks and derive the fraction of the observed sample that will escape the Galaxy. We discuss 
the form of the velocity distribution function and give arguments that suggest that kicks 
>50 km s~ x always occur and that a significant fraction of objects are moving so fast that 
only strong supernova kicks can explain their motion. We speculate about the relationship 
of the shape of the velocity distribution function to binary evolution scenarios. We discuss 
the implications for core collapse, drawing attention to the breaking of symmetry implied 
by the observed kick size. 



Finally, we include a discussion of the search for the highest velocity objects in §11 



We summarize our results in §12. Individual objects with idiosyncratic distance estimates 



or with extreme birth velocities or z altitudes are discussed in Appendix 



2. PULSAR SAMPLE 

The observables we use to analyze pulsar space velocities are the spin period and its 
first time derivative (P and P); galactic coordinates (£,b); the dispersion measure (DM); 
and the proper motions in the directions of right ascension and declination (/zra, /^dec) an d 
their associated measurement errors (<7 MJL4 , cr^ DBC ). We incorporate auxiliary information 
on particular lines of sight that includes the scattering measure (SM) and the locations of 
HII regions that may perturb distance estimates based on DM. The distance scale is based 
on the electron density model of TC, augmented by neutral hydrogen (HI) absorption 
measurements and other independent constraints on the distances, including 2 parallax 
measurements. In aggregate, the available data provide a formal distance estimate and, for 
most objects, lower and upper bounds (D L) Du) such as those presented by Frail (1990), 
Frail & Weisberg (1990) and Koribalski et al. (1995) based on HI absorption measurements, 
pulsar associations with supernova remnants, and parallax measurements. 

Table [l] lists the data that we have used, which results from the requirement that the 
spindown age (P/2P for braking index n = 3 and no torque decay) be less than 10 Myr and 
that the distance estimate include a lower and upper bound. A chronological age cutoff is 
necessary because we treat the orbital motion of the pulsar analytically by a Taylor series 
approximation and this ceases to be accurate once significant changes in the acceleration 
have occurred. We focus on the vertical pulsar motion in our analysis ignoring all variations 
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with Galactocentric radius. We choose the age cutoff of 10 Myr to guarantee that virtually 
all bound objects are still in the first quarter cycle of their oscillation perpendicular to the 
galactic plane. As we discuss in more detail below, the 10 Myr age cutoff assures quite 
accurate vertical solutions in an infinite, planar Galaxy. However, in 10 Myr a 1000 km/s 
pulsar moves ~ 10 kpc, so the assumption of uniformity in the plane is nontrivial. 

For each object the chronological age is based on P, P and a spindown model. It turns 
out that the models of greatest interest have n > 3 and/or torque decay. For them, the 
chronological age is less than P/2P so the imposed cutoff of 10 Myr is more conservative 
than need be the case. Even for 2 < n < 3, the maximum chronological age of 20 Myr 
degrades but does not destroy the accuracy of the vertical solution. The age cut implies 
that our sample excludes millisecond pulsars and other pulsars that may have undergone 
accretion-driven spinup, a process that necessarily is slow and is therefore manifested in 
older objects. 

Of 96 objects with proper motion measurements, 47 satisfy our age and distance 
criteria. The 49 excluded objects include 45 that are older than 10 Myr, of which 18 are 
older than 100 Myr. An additional 4 objects were excluded because the TC distance model 
yielded only a lower bound on the distance and there was no auxiliary information to 
provide better constraints on the distance. The sample of 47 objects includes three with 
parallax measurements (B0630+17 [Geminga], Caraveo et al. 1996; B0823+26, Gwinn et 
al. 1986; and B2021+51, Campbell et al. 1996). 

3. RADIAL VELOCITIES AND BIRTH ALTITUDES 

We first derive constraints on the radial velocity and birth altitude, zo, from 
measurements of proper motion and distances on individual objects. In so doing, we 
demonstrate that all pulsars with measured proper motions are consistent with birth near 
the galactic plane, \z \ < 0.3 kpc. This kind of analysis was first applied by Helfand & 
Tademaru (1977) to a sample of 12 objects. Our conclusion about pulsar birth in the 
galactic disk is identical to that of Helfand & Tademaru. Similarly, we find that the 
apparent motion of a few pulsars toward the galactic plane is due exclusively to projection 
effects. We extend previous analyses by taking into account errors on estimated distances 
and proper motions, acceleration in the galactic potential and alternative values for the 
braking index. 

In the following, we review the age uncertainty and describe a model which incorporates 
an arbitrary braking index and torque decay. We outline the geometric relations that exist 
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between the estimated distance, proper-motion components, z-distance from the Galactic 
plane, z-value at birth, radial velocity and pulsar age. We formulate a likelihood function 
for a quantity composed of possible z-values at birth and peculiar radial velocities today. 
We compare its value on a case-by-case basis to a priori expectations. 



3.1. Age Estimate Uncertainties 

The spindown age for an object with period P and period derivative P is defined as 

Ts " (2) 

(n — 1)P 

where n is the braking index. We have n = 3 for braking by radiation from a point magnetic 
dipole of constant strength. If n is a known constant, n > 1, and the initial period P <^ P 
then the spindown age t s and the chronological age t agree. 

Measured braking indices are smaller than the dipole result: n = 2.28 for PSR 
B0540-69 (Boyd et al. 1995); n = 2.51 for the Crab pulsar (Lyne, Pritchard & Smith 1988); 
and n = 2.83 for PSR B1509-58 (Manchester et al. 1985; Kaspi et al. 1994). Most pulsars 
are much older than these objects and their braking indices cannot be measured reliably 
because timing irregularities mask the contribution from the small expected second period 
derivative, P. Work below indicates that values of the braking index n ^ 3 (and no torque 
decay) are less consistent with the kinematic data than are values ^3 (and/or torque 
decays on a time scale of a few million years). Since most of the pulsars in our sample are 
significantly older than those for which braking indices have been measured directly, these 
facts imply that braking indices are not constant and/or torque decay is significant over a 
pulsar's lifetime. 

In principle, many physical effects can alter spindown from the static magnetic dipole 
limit. The strength of the magnetic field may change: Ohmic decay (Ostriker and Gunn 
1969) and field growth (e.g. by thermal effects, Blandford, Applegate and Hernquist 1983) 
are both possible. Particle emission processes may be more important than dipole radiation 
in angular momentum loss (Michel 1969, Goldreich and Julian 1970) and the losses might 
scale differently with frequency and with magnetic field than dipole losses. The dipole 
geometry may change in time: some systematic changes in the relative orientation of the 
spin axis and the dipole axis were observationally inferred from the shape of pulsar beams 
(Lyne & Manchester 1988). The effective size of the dipole may change in time (Melatos 
1997). Surface fields (as opposed to the field at the light cylinder) may be altered by mass 
accretion (Romani 1990) and crustal plate motion (Ruderman 1991); whether such changes 
influence the spindown rate is unclear. 
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It is difficult to test these ideas directly since measurement of torque decay may be 
masked by other physical processes occurring inside the neutron star. Some information 
is available. General observational indications are that pulsars that are or were members 
of binaries have weak fields of 10 8 — 10 9 G that do not decay over lifetimes of 10 9 — 10 10 
yrs. In a number of recent, specific analyses involving accretion in binary systems field 
decay has been inferred to occur (e.g. Burderi, King & Wynn 1996). Current statistical 
studies of pulsar properties do not find strong evidence for field decay in single pulsars 
(e.g. van den Heuvel 1993) nor for changes in the relative orientation of spin and dipole 
axes (Bhattacharya & van den Heuvel 1991). Recently, Wang (1997) has argued that an 
isolated pulsed X-ray source accreting from the interstellar medium must have undergone 
field decay from ~ 10 12 G to 10 10 G in ~ 10 7 yrs (by power law decay) or ~ 10 s yrs (by 
exponential decay). 

In view of the uncertainties, we adopt a phenomenological approach and write 
Q = —K(t)Q n , where t is the chronological age and K(t) describes any time dependence of 
the torque law other than a pure power of the spin rate ("torque decay"). Most discussions 
on torque decay refer explicitly to exponential decay of the magnetic field for n — 3. Then 
K(t) oc B 2 , and the quoted field decay time is twice that of the torque decay. The notion of 
magnetic field decay has changed considerably since it was first introduced in discussions 
of pulsar evolution. Ohmic decay in the neutron star crust was initially thought to be fast 
enough to account for the turn-off of radio emission from pulsars (Gunn & Ostriker 1970) 
but it is now regarded as too slow to be relevant (e.g. Goldreich & Reisenegger 1992). For 
our uses, the phenomenological model applies for any n > 1 but the relation between B 
field and K(t) will vary. 



For exponential torque decay, K(t) oc exp(— t/r K ) and the chronological age is given by 
(Bailes 1989) 



If t <C tk and P ^> P , then t ~ rg. In other words, the spindown age and the chronological 
age agree when the torque decay has not taken place and when the initial period is short 
compared to the current period. On the other hand, if t > tk or if P ~ P , then t <^ t$ 
implying that the spindown age overestimates the chronological age. 

A braking index may be inferred from the instantaneous rate of change of frequency. 
The estimated braking index is 




(3) 




(4) 



For Ts <C tk, n = n but for ts <; tk and n > 1, n > n. This sort of evolution would 
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be consistent with both the small measured values of n for very young pulsars and the 
larger kinematically inferred values of n for older objects (see below, § |4.5| , where we discuss 
constraints on spindown laws). 



3.2. Coordinates, Kinematic and Peculiar Velocities, & Proper Motions 

We adopt U, V } W coordinates (Mihalas & Binney 1981, pp. 382-383) where U and V 
point toward i = 180° and I = 90°, respectively, for b = and W points toward b = 90°. 
We transform the measured proper motions and their errors to the UVW system to obtain 
^u,v,w and o mvyif . 

Let z be the present distance from the plane, zq be the distance at birth from the 
plane?] and z be the acceleration at the current position. The present-day, "kinematic" z 
velocity of an object of age t is 

t/ Z — Zq Zt 

V z (5) 

t 2 K ' 



and the birth z velocity is 



V z0 &V Z - zt. (6) 



The Taylor series expansion is accurate for pulsars younger than ~ 1/4 their z-oscillation 
period. Figure [l] shows V z (t) for pulsars with zq = moving in a three-component model 
for the galactic potential (Paczynski 1990). The very slowest pulsars (e.g. V z0 < 10 km 
s _1 ), have asymptotic periods of 60-70 Myr (near the Sun), while somewhat faster pulsars 
(V z q ~ 50 km s" 1 ) have periods ~ 100 Myr. By restricting our analysis to objects with 
Ts < 10 Myr, we can safely assume that our expressions above accurately account for 
vertical acceleration. For example, with V zQ = 100 km s -1 at 20 Myr, the Taylor series 
calculated velocity differs from the true current velocity by less than 3% despite changes 
of order 33%. Faster pulsars show smaller relative errors and slower ones, larger errors. 
(For initial z velocities of 20 km s _1 , the error is 4 km s _1 in the final z velocity.) Since 
our main purpose is to describe pulsar velocities statistically over the broad range of values 
encompassed by pulsars, this precision is more than sufficient. As a practical matter, 
errors in pulsar distances and the proper motions themselves dominate uncertainties in the 
approximations associated with the vertical motion. 

The pulsar peculiar velocity (motion with respect to the mean disk in the vicinity of 



2 We assume the Sun is at midplane, even though there is evidence that it currently is ~ 15 pc above the 
plane (Hammersky et al. 1995). 
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the pulsar) is 

V( p ) = Dfj,^ + K (P) fi (7) 

where D is the distance, is the peculiar proper motion and V r ^n is the line-of-sight 
peculiar radial velocity. The proper motion = + yS DGR ) where ^ DGIV > is the known 
contribution from differential galactic rotation (DGR). The actual radial velocity is also 
composed of a peculiar part and a known contribution from DGR (V r = V,P + V r ^ DGR ^). 
The z-component of the peculiar velocity is 

sin b (8) 

and V z = VP since the effects of rotation lie in the plane. One combination of the 
individually unmeasurable quantities V( p > and zq is determined in terms of the distance, 
current acceleration, peculiar proper motion and age (assumed known): 

t sin b t sin b V 2 J 
Note that V r e ff relates a particular combination of initial birth location and current 
peculiar radial velocity to directly measured quantities. 



3.3. Likelihood Function for Radial Velocities & Birth Altitudes 

We will use K, e // to characterize individual objects: large values imply large peculiar 
radial motions and/or large initial birth heights. We seek the pdf of V T ,eff given the 
uncertainties in the proper motion measurements and in the distance to individual objects. 
We treat this analysis as a preliminary attack on the determination of the 3D velocity 
distribution which is the subject of the following section. 

We write the probability density function for the error in the proper motion 
measurement fsu,(S/J,) where the error is 

5/j> = /x — /2 (10) 

and fx is the reported value. We assume that the pdf for 8/a is a Gaussian with 
reported standard deviation a^. For the distribution of K,e// we are concerned with the 
component of proper motion out of the plane. Denote by K-,e// the value of K,e// when 
l^w = fiw — fJ-w ^ hi Eq- H above. The pdf for an object with distance D, in direction n 
and chronological age t is 

fv r ,eff( v r,eff I (j>w,D,n,t) = J d5n w f Sflw (5fi w )5(V rteff - V r , e ff - [D5/j, w / sinfe]) 

sin b (p) (dgr) ~ \ 
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(p) 

In the second line, [i w i s given explicitly by 

(p) sin b 'zt sinfo , . 

M = — + 2D- ~D~ (12) 

where Hw GR ^ and z are the position dependent differential galactic rotation contribution 
to the proper motion and the position dependent galactic acceleration. (For notational 
simplification, we suppress writing the direction and age in the pdf's that follow in this 
section.) 

There is significant uncertainty in the distance to most pulsars. The pdf of K,e// is 
expressed in terms of the range and distribution of distances allowed by observations. We 
have 

fv r , eff (V r , eff \fi w ,D*) = JdDf D (D)f Vreff (V r>eff \» w ,D) (13) 

where D* stands for the imperfect distance knowledge and fo{D) is the corresponding 
pdf for distance. The specific form for the distance distribution may not be critical for all 
applications. If we are primarily looking for objects that are grossly discrepant with what 
we know about the disk scale height and the typical velocity of pulsars today, then a simple 
approximation may be sufficient. We assume a flat distribution for the distance in the 
interval [Dl, Djj] where the lower and upper distance bounds are based on considerations of 
the TC distance model along with HI absorption and other measurements. In later sections, 
we will use an expression for Jd{D) that is based on the expected number of pulsars in a 
population that would be visible. 

The pdf has the form 
fv r , eff (V r ,eff I frw,D*) = Sm& / dD D~ x fs^iip +Hw° R) -M- (14) 

We use it to derive estimates for the covariant combination of the peculiar radial velocity 
and birth altitude of a given pulsar. The pdf determines the most likely value for V rie ff and 
associated confidence interval \V r ,eff~ K,e// + ]- Then values for the peculiar radial velocity 
V^ p ) and birth height Zq and individual confidence intervals may be found by using Eq. |9|. 

Initially, our analysis uses the chronological age t = t$ for braking index n = 3 and 
no torque decay. Later in the paper, we consider alternative age estimates. To calculate 
the accelerations 'z in Eq. |l^ we use a three-component model for the galactic potential 
(Paczynski 1990). To model differential galactic rotation, we use a flat rotation curve with 
circular velocity v rot = 220 km s -1 . Figure || shows plotted against —0.5 < z < 0.5 
kpc for all the pulsars in the sample. Also shown for three pulsars are dashed, parallel lines 
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that designate the 68% likelihood range in due to measurement errors in proper motion 
and to distance uncertainties. The figure demonstrates that, for most objects, solutions for 
(y( p \ Zq) exist with z < 0.3 kpc and |V^ p )| < 10 3 km s _1 . Three objects (the Crab pulsar, 
B0540+23 & B0611+22) require \z \ ~ 0.2 - 0.3 kpc within the plotted range of V r {p) . The 
allowed solutions imply that (1) the pulsar data are consistent with birth in the plane from 
Population I progenitor stars; and (2) no known pulsar younger than 10 Myr with a proper 
motion measurement need be born high above the galactic plane. 

We also considered alternative relations between spindown and chronological ages. 
The braking index value n = 2.5 yields a roughly similar set of solutions but the situation 
becomes increasingly extreme as n decreases. For n = 2 some objects (e.g. 1933+16) have 
large K. j6 //, e.g. |V r (p) | > 5000 km s _1 at z = kpc, or \z \ > 0.3 kpc at V^ p) = km s -1 . 
Figure |3] shows (V} p \ zq) solutions for four pulsars for braking indices of 2, 2.5, 3 & 4. 

4. THE 3D VELOCITY DISTRIBUTION 

By combining the radial velocity constraints of §[5] with proper motion measurements, it 
is possible to derive information on the full 3D velocities of pulsars and their birth distances 
from the galactic plane. Here we derive the joint pdf for birth z and birth velocity. 

We first relate the initial distribution function for pulsars at birth to the observed 
distribution function today. Then, we derive a likelihood function for proper motion data 
to be used to infer information about the initial peculiar velocity and birth height above 
the plane. 

4.1. Birth & Present-day Distribution Functions 

The distribution function for particles moving in a fixed background potential is related 
directly to the distribution function of initial conditions. Choose an inertial frame. At time 
t let the position be X, the velocity V and let the corresponding initial conditions be X 
and Vo- Liouville's theorem states that the final and initial distribution functions satisfy 
/x,v(X, V) = ixo,v (Xo, Vo). 

We will apply this result to relate the initial and final pulsar distributions, bearing in 
mind that V is the sum of the circular rotation velocity and the peculiar motion, including 
the birth kick. The distribution of initial conditions /x ,v (Xo, Vo) is expressed in terms of 
the initial peculiar velocity distribution / (p) (X , V P ^). If its age is sufficiently small, 
then the equations of motion show that the changes in peculiar velocity are small. To be 
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more precise, we ignore terms of 0(Qt, Vq t/R) compared to V /v ro t where Q is the 
circular frequency, R is the Galactocentric radius and v rot is the rotation velocity. We 
also ignore variations with R, so we factor /x ,v = fgp x fz,v ( z -> Vq) into a constant 
two-dimensional, galactic plane part (f gp ) and a four-dimensional pdf in z and Vo. This is 
justified for particles (neutron stars) whose total radial displacement is small compared to 
the length scale for variations of disk properties. This approximation may be problematic 
given the maximum distances traveled in 10 Myr. The physical content of these assumptions 
is that the birth rate, birth scale height and initial velocity kick distribution do not vary 
appreciably over the range of the radial motion, at least for parts of the Galaxy that 
contribute to the observed pulsar sample. We find from Eqs. [7| and |8| that 

z = z - V z t + X -zt 2 = z- {Duff + V r {p) sin b)t + ht 2 (15) 

v£ p) = V( p ) - zzt = (ZV P > + V} p) n) - zzt. (16) 
The pdf of present-day pulsars is related to the pdf at birth by 

/x,v(X, V) = / X0)V (P)(X , VS P) ) = / S p4 iV (p)(%vl P) ), (17) 

For convenience, we abbreviate / (zq, Vq P ^) = f(zo, V^). 



4.2. PDF of Proper Motion 

Now we derive the likelihood function for the proper motion. We begin by forming 
the likelihood of a perfectly measured proper motion /j, for a pulsar with perfectly known 
position X, displaced from the observer's position by Dh. The pdf is 

///MX) = j^j rf 3 v/x,v(x, v)5V - (V ~J o)± ). (18) 

Here / x ,v(X,V) is the present-day pulsar pdf; A/Vt|x = / ^ 3 V/x,v(X, V) = /x(X) 
normalizes the proper-motion pdf; and V is the velocity of the Sun due to galactic 
rotation. 

The distribution of the measured proper motion jl, including measurement errors, is 

/ A (A|X) = J d/*/(£|Ai)/Ax(/i|X), (19) 

where f(fl\fi) is the likelihood of obtaining the measured value given the true one and is 
simply the usual Gaussian function for each component with standard deviations that reflect 
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measurement uncertainties. From our previous discussion, we assume f(ft\n) depends only 
on the difference, Sfi = (J, — /2, written as fsu,(S/J,). Let D* stand for lower and upper cutoffs 
for D based on dispersion measure, known HII regions and so forth. Let fu{D\D*, n) be 
the resultant probability distribution. Marginalizing over the true distance, 

/ A (£|£>*, n) = J dD f ^(fl\X.) f D (D\D* , n). (20) 

We can write the pdf for distances in the pulsar population described by /x,v assuming all 
such objects are visible. We find 

M D\D>,t) = J £g£L ry (21) 

where the integral is limited to the range implied by D*. This is the point where we 
have ignored the detection-related aspects of the pulsar problem. Integrating over the 
perpendicular velocity we find 



/ A (£|L>*,n) = Wrf^ J dDD 4 J d» J dV r f s ^ ~ aO/x,v(X, V) 



(22) 

V=V +D^+y r n 



where V r is the radial velocity in the inertial frame and Nrf D * fi normalizes the pdf. To 
carry out the integrations (dfj,dV r ) we change variables as follows: — fj,— /x^ DG ^(n, D) 
and =V r - Vy GR \n, D) and use Eqs. [15|-[17|to give 

f fi fr\D;&,t)=tfrf D . A J dDD 4 J d^ J dVyf 5fl (»-» {P) -» {DGR) )f(zoM P) )- (23) 

Note that and are dummy integration variables and that the transformation to 
the initial peculiar velocity distribution requires knowledge of the chronological age. The 
data and the DGR contribution to the proper motion occur only in the pdf that describes 
the error distribution in the form of the difference: ji — ^ DGR ^ . For cases of interest (ts ^ 10 
Myr), the acceleration correction to the velocity argument in Eq.[l6| is only a few tens of km 
s _1 and is, therefore, unimportant for establishing the velocity distribution of high-velocity 
objects. However, the correction to the spatial argument in Eq. [15] is significant. In the 
following, we ignore the velocity correction but apply the spatial correction. 



4.3. Parameterization &c Likelihood Function 

The likelihood function for the parameters is the product over N psr pulsars, 

Npsr 

C=U£k, (24) 
fe=i 
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where the likelihood factor, £ fc = fn(fi k \xik, D%, t^), for each pulsar is simply Eq. |23 
evaluated using the measured proper motions (and errors) along with the distance 
constraints D L , Djj, the direction n and the age t. 

To apply Eq. ^3], we adopt a parametric approach, where we assume particular forms 

(p) (p) 
for the pdfs in z an d V . Specifically, we assume that Zq and V have a multicomponent 

Gaussian pdf of the form 

n g 

f(zo, Vj p) ) = J2 Wjgid(z , h zj )g 3 d(v£ P \ a Vj ). (25) 
i=i 

In Eq. ^5| gidip, cr) is a standard ID Gaussian pdf with zero mean and standard 
deviation a, gi d {p,a) = (27T(t 2 ) _1 / 2 exp(— p 2 /2a 2 ), while g 3 d is a 3D Gaussian 
g3d(q, cr) = (27rcr 2 ) _3 / 2 exp(— q 2 /2a 2 ). The weights Wj sum to unity. The pdf is 
defined so that fdzod^V^ is the infinitesimal probability. 

We have chosen a multicomponent Gaussian model because its analytical properties 
allow it to fit a wide range of shapes for the actual distributions in zq and V^ p '. There 
is not necessarily an implied physical basis for this choice of form: the different Gaussian 
velocity components need not correspond to different population components. 



4.4. Results & Comparison of Models Using Odds Ratios 

The parameters to be determined for a pulsar sample are (1) n g standard deviations for 
velocities, avj', (2) rih < n g scale heights for the birth altitude, h z j; and (3) n g — 1 weights, 
Wj for a total of n g + rih — 1 parameters. 

We considered a set of models with increasingly complex velocity and birth height 
distributions. We label the models by n g .rih and briefly describe them below: 

1. [Model 1.1]: a single component model (n g = = 1) 

2. [Model 2.1]: a two component velocity model with a single scale height (n g = 2, nh = 1) 

3. [Model 2.2]: a model with two scale height and two velocity components 
{n g = 2, rih = 2) and 

4. [Model 3.1]: a three component velocity model with a single scale height 
(n g = 3,n h = 1) 
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We assume flat priors for the parameters in selected ranges that are listed in Table 0. 
We generated the posterior probability distribution of parameters which is just the 
evaluation of the likelihood in the selected ranges. The modes of the posterior, i.e. the 
maximum likelihood results, and log £ appear in Table [| These are the "best" values of 
parameters for each model. The analysis was repeated for several age cuts for the sample: 
T Smax = 1, 5 & 10 Myr (for n = 3 and no torque decay). The mode of the distribution is 
typically well-defined. For example, in Figure f| we display contours of log C (log likelihood) 
for Model 1.1 with Ts max = 10 Myr. The mode is identified by the cross. In Figure || we 
display contours for Model 2.1 in two-dimensional surfaces that pass through the mode. 

The results in Table [3] clearly suggest the presence of a significant high velocity 
component when the models are sufficiently complex. At the same time, examination of 
the likelihood contours in the simplest model does not provide any obvious indications that 
its single 300 km s _1 component is inadequate. To make such a determination, we must 
compare the models in a systematic manner. 

To compare models we use the odds ratio (cf. Gregory & Loredo 1992) to quantify 
goodness of fit while penalizing models with more parameters. We assume that all models 
are equally probable, a priori. The odds ratio reduces to the "Bayes factor," which is 
the ratio of global likelihoods of two models (Gregory & Loredo's Eq. 2.12). The global 
likelihood for a model M given data T> is 

f(V\M) = J dOf (O\M)£(O), (26) 

where fg(0\M) is the prior pdf for model parameters 6 and C(0) is the likelihood function 
as we have used it in this paper. We have already assumed the prior pdf is flat and we have 
seen that C{0) is strongly peaked. Thus, letting 6 be the parameters that maximize C(0), 



f(V\M) ~ f e (0\M) f dOC{0) = ^, (27) 



A/ 

where the last equation defines the integrated likelihood, A M = J d0C(0), and the volume 
in parameter space that is searched, Vm = fa(0\M) 

We take the single Gaussian model with two parameters as our reference model: 
Mi = 1.1. The odds ratio for the M th model relative to M\,\ becomes 

f(V\M) _ A M V Ml 

which we evaluate through numerical integration of the likelihood function for each model 
over a uniform grid in parameter space with bounds given in Table 0. The associated 
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volumes for the four models are: 



Vi.i = Ah zl Aa Vl 

V 2 .i = A^iA(Ty2^2.1 
V3.I = V1.1 Aw7iAw 2 A(Jy2Acrv3v43.i 



(29) 



where ^2.1,^2.2 and A3.1 < 1 are factors that account for the overlap in our search of 
parameter space for the maximum likelihood and the resulting degeneracy of the models. 

The maximum likelihood values and the volumes yield the odds ratios given in Table [| 
The odds ratios show that the data strongly favor the multiple-component models (models 
2.1, 2.2 and 3.1) over the single-component model (1.1). However, the additional complexity 
of models 2.2 and 3.1 over 2.1 is not demanded by the data. This conclusion holds for the 
two largest age cutoffs considered (5 and 10 Myr). For a 1 Myr cutoff, the small number of 
data points appears to reduce the discriminating power between models and, consequently, 
we have ignored these results. The model fitting suggests that a single birth-z scale height 
of 0.13 kpc applies to both velocity components, which have rms velocities (3D) ~ 175 km 
s _1 and ~ 700 km s -1 . 

To refine the estimation, we derive the marginal distribution of each parameter. 
For a given parameter 9j G 6, the marginal pdf is the normalized integral over all other 
parameters 

L xce .dO£{0) 

f *M= TkW) ■ (30) 

where the integral subscript 'exc. 6/ means that all parameters except the ] th one are 
integrated over. Fig ^| illustrates the marginal distributions for Model 2.1. Table f| gives 
the best-fit values of parameters for the preferred Model 2.1 along with 68% confidence 
intervals. The confidence interval was calculated by finding the region in each marginalized 
distribution containing 68% of the area. 



^jWjUvj ~ 264 km s" 1 and the three mean speeds are 211, 331 and 421 km s" 1 . The 
speeds for the three-component Gaussian model (3.1) are nearly identical, while they are 
about 14% larger for the single-Gaussian model (1.1). 




The mean ID, 2D and 3D speeds are given by (Vm) 
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4.5. Spindown Ages & the Braking Index 

The likelihood analysis we have presented is based on a braking index n = 3, an infinite 
torque decay time, t k (cf. Eq. |3|), and a negligible birth period, P . The analysis is sensitive 
to the spindown model because we infer the chronological age t from the spindown age ts- 
For example, Eq. |23] shows that the likelihood function depends on the initial distribution 
at z Q and Eq. [H| in turn, shows that the birth position zq is a direct function of t as well as 
other observed (z) and marginalized {V^ p \ ityp) quantities. Here, we investigate how the 
results for the 3D velocity distribution depend on the braking index and decay time. As 
mentioned earlier, measured braking indices for very young pulsars (age ~ 10 3 to 10 4 yr) 
are less than 3, yielding spindown ages that are larger than we have assumed in our use of 
the conventional spindown age with n — 3. As n and tk are varied, ages for some objects 
increase while others decrease. 

Figure [7] shows the log likelihood for the double Gaussian model (2.1) plotted against 
torque decay time, tk, for various values of braking index. All model parameters were 
varied to find the maximum likelihood. The likelihood peaks at (n, tk) ~ (2.5, 3 Myr) and 
the range of variation in the final likelihood value is small compared to that due to other 
model parameters. For example, the likelihood at (n, tk) = (3, 5 Myr) is only 22% smaller 
than the peak. If torque decay does not occur (i.e. — ► oo), the peak likelihood is found at 
n rs 4.5, but with a likelihood that is a factor 20.3 smaller than for (n, tk) ~ (2.5, 3 Myr). 
If torque decay is due to magnetic field decay with n = 3, the field decay time is twice the 
torque decay time, or ~ 10 Myr. 

The best values for the birth scale height, velocity parameters, and relative weights 
do not change appreciably with braking index and decay time. (It should be pointed out 
that we have used a fairly coarse grid in this analysis: Aw\ = 0.1, Ah zl = 0.05 kpc, 
Aovi = 50 km s _1 and Aoy 2 = 200 km s _1 .) On the other hand, the mean sample 
age does vary in a easily characterized way. It is 2.95 Myr for (ti,Tk) = (3, oo). For the 
best fit with (n, tk) = (2.5, 3 Myr) it is 1.64 Myr while the best fit with no torque decay 
for n = 4.5 yields 1.68 Myr. In the case with no torque decay, the best fit reduces all 
age estimates by the same factor, 3/4.5 = 0.67. By contrast, for the absolute best fit for 
(n,r K ) = (2.5,3 Myr), some age estimates are increased while others are decreased. 

Even though all model parameters were varied, the best parameter values for 
wi, h zl , ay i and oy 2 do not change significantly as the braking index and torque decay time 
are varied. Certain model parameters (e.g. oy) are directly tied to a measured quantity 
(n^) and, hence, not susceptible to changes in t. Other model parameters (e.g. h z ) are 
strongly influenced by young objects at small sin b that have not had the opportunity to 
move far for any plausible ages. In addition, there are substantial errors in distance and 
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proper motion measurements and these overlay the intrinsic timescale change. Typical 
distance errors are factors of ~ 2 and typical proper motion errors are ~ 50%. Thus it is 
not unreasonable that we distinguish spindown models that have typical age estimates that 
differ by ~ 2. 

The pulsar sample we have analyzed is more consistent with ages that are smaller, 
on average, than the conventional spindown ages. If no torque decay occurs, then "large" 
braking indices, 3<n<6 are preferred over small braking indices, such as n ~ 2.5 measured 
for the Crab pulsar and n ~ 2.4 for B0540-69 in the LMC. If torque decay occurs, as is 
suggested by our best fit model, then a braking index of 2.5 — 3 is preferred. An average 
braking index as small as 2 yields a likelihood that is ~ 10 2,9 times worse than the best fit 
model. 

Our results are therefore at odds with those of Lyne et al. (1996), who estimate a 
braking index n ~ 1.4 ± 0.2 for the Vela pulsar and who speculate that braking indices 
decrease as they age, at least in going from Crab-like pulsars ~ 1000 yr old to Vela-type 
pulsars about 20 times older. Lyne et al. conclude that pulsars are generally older than the 
conventional spindown age, in conflict with our conclusion. 

Assuming the braking index for the youngest objects has been accurately determined, 
our results could be consistent with one of several proposed scenarios in which the braking 
index increased with age. Finite size corrections to the dipole spindown model (Melatos 
1997) predict such an evolution. Heintzmann & Schrueffer (1982) noted that braking 
indices of older pulsars should be 3 or more based on plasma effects in the magnetosphere. 
Muslimov & Page (1996) suggested that, if the crustal magnetic field of neutron stars grows 
at early times (0.1 - 10 kyr), then braking indices should increase from n < 1 and asymptote 
to n — > 3 on a time scale of ~ 10 kyr. 

In this analysis we have assumed that the birth spin period is much smaller than the 
present day period. There is evidence of discrepancies between spindown ages of some 
millisecond pulsars and ages constrained by temperatures of their white-dwarf companions 
(e.g. Camilo et al. 1994). Vivekenand & Narayan (1981) and Phinney (1996, private 
communication) have argued that pulsars may be born with nonnegligible spin periods 
("injection"), in which case spindown ages overestimate chronological ages. The fact that 
the mean sample age of our most likely models is less than that of the magnetic dipole 
model without torque decay would be consistent with this suggestion. However, it does 
not follow that our analysis favors injection. In varying the spindown laws, the ages of 
all objects change, whereas if pulsars have a range of birth spin periods only objects with 
current periods comparable to birth periods have significantly altered ages. The set of 
spindown models we have tested does not cover the hypothesis of injection. 
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In we explore the issue of pulsar ages in generality by deriving the distribution of 
chronological age for individual pulsars. This treatment starts from the kinematic data 
for each pulsar and our inferred scale height and velocity dispersion to deduce the age. 
The treatment is independent of the spindown model, in so far as the scale height and 
velocity dispersion parameters are themselves insensitive to the precise form of the spindown 
model, as appears to be the case. The general conclusion of the case-by-case analysis is 
consistent with the above conclusions: pulsars are on average roughly 50% younger than 
their conventional spindown ages. 



5. 3D VELOCITIES OF INDIVIDUAL OBJECTS 

In this section we use the global properties of the observed pulsar population (from 
§3) to aid estimation of the properties of individual pulsars. First, we write out a formal 
expression for the pdf for V^ p ) and z$ for an individual object given its distance limits, 
proper motion measurements, and the chronological age. In constructing the pdf, we 
incorporate as prior information the previously derived parameters for the distribution 
of birth altitude and velocity. Second, we Monte Carlo the resulting pdf to determine 
mean and mean square values for each pulsar in the sample. This highlights a number of 
interesting objects. 

The analysis here is similar to the analysis of §|3] but differs in that we estimate 
transverse speeds for each object and incorporate what we have learned from the 3D 
analysis of the whole population. The 3D velocities of individual objects can be used to 
test for correlations with other pulsar parameters. 

Let Q be a quantity whose pdf we seek; Q may be a function of any of the observables, 
with functional form Q. Signify the data for a given pulsar as T> = {D*, /2, n}. For each 
object we estimate the chronological age t from P, P and an assumed spindown model. Let 
I represent the background information, including the posterior probability distribution for 
the model parameters governing the birth scale height, the velocity kick distribution and 
the relative weights of components. Since the posterior was found to be strongly peaked at 
the mode and the mode insensitive to the spindown model, we simply assume that X sets 
the model parameters at the mode (i.e. the most likely values). We derive fq(Q\V,t,T) 
following the same general procedure used to derive Eq. First, form the pdf for Q in 
terms of the distribution function for today's pulsar positions and velocities, assuming no 
observational errors in fi and distance. Second, marginalize over the distance range and 



the proper motion errors. As before, we take Eq. 21 for the pdf for distance D between 



the upper and lower limits. Third, make the change of variables d 3 \ — > D 2 d^S F " > dV^ p " > 
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and replace the present day pulsar distribution function with the corresponding initial 
distribution function. The result is 

f Q (Q\V,t,X) dDD'J V P) / dVrfsnfr-nW -^ DG ^)f{z Q yP)5{Q-Q) 

(31) 

where Mq 1 is the normalization constant; z$ and Vq^ are the initial conditions (given 
by Eqs. [15] and [l6j) ; and the range of D integration is determined by D*. The previously 
derived parameter values for birth scale height and peculiar velocity and the age occur in 

As a specific example of immediate interest, we take Q = {zo,V^} to give the joint 
pdf for the present day peculiar velocity and the birth altitude of an individual pulsar 



f zo ,v<Mzo,V {P) ) 



oc 



dD 

dZn 



D = D 



z + z • V( p ) - \zP , 

D = — 2 - . 32 

sin b 

The solution for D is implicit because the vertical acceleration depends on D. 

The expression clearly shows the role of the initial birth and velocity distribution 
f(zo, Vq ) on limiting the allowed range of properties for a pulsar observed today. For 
poor or no observational constraints on proper motion, the birth distribution functions as 
a prior. However, if a proper motion measurement is precise, then fsfj, provides a strong, 
perhaps dominant constraint on the velocity components and birth z for the object. 

To apply Eq. |32| for each pulsar, we generate 200 Monte Carlo samples of Zq and V( p ) 
using the method of rejection (cf. Rubinstein 1981; Press et al. 1992). For the pulsar age, 
we have assumed the canonical age given by Eq. with n = 3 (i.e. magnetic dipole radiation 

and no torque decay). Table |5| gives the mean and rms values for V = \j V r {P)2 + \V { p \ 2 , 

V^ p \ V^, and z for the full sample. The two components for the perpendicular 
velocity (V±\, Vj.2) given in the table correspond to components in the right- ascension and 
declination directions. The rms values indicate the spread in possible solutions for each 
variable, subject to measurement errors in the proper motion, distance estimation errors 
and the fact that only two velocity components are determined by proper motion. Most 
objects show small average values for V}- p \ but with a large spread due to uncertainties 
in distance and proper motion. On the other hand, three objects require significant radial 
velocities, V r ^ > 3a v ( P) : B0943+10, B1953+50 and B2154+40. The latter two objects 
require radial velocity magnitudes ~ 1500 and ~ 1000 km s -1 , respectively. There are 
8 objects with 3D mean velocity estimates exceeding 500 km s _1 : B1508+55, B1706-16, 
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B1842+14, B1933+16, B1946+35, B1953+50, B2154+40 and B2224+65. Of these, five 
have velocities in excess of 1000 km s _1 . Finally, 9 objects have birth z values that are 
determined to better than 3a Zo : B0136+57, B0531+21 (Crab), B0540+23, B0611+22, 
B0630+17 (Geminga), B0656+14, B0740-28, B0833-45 (Vela), and B1929+10. Of these, 
three pulsars have |zo| ~ 0.2 kpc. It is notable that these high z pulsars are in the galactic 
anticenter. 



6. PULSAR AGES 

Now we derive the pdf for the chronological age, t, of a pulsar given its proper motion 
and distance measurements. Previously, we estimated t from the spindown time ts or its 
variants using alternative braking indices and torque decay times. Here, we rely solely on 
the kinematic information associated with each object and with the global parameters that 
characterize the birth scale height and initial peculiar velocity. These parameters were 
determined for a range of spindown models in § [4.5| and the results were found to be only 
weakly dependent on breaking index and torque decay time. This simplification means that 
the pdf for t is insensitive to the spindown assumptions made in the inferring the scale 
height and peculiar velocity parameters. 

As we have previously emphasized, uncertainties in position and velocity are 
inextricably linked to the age uncertainty. We use Bayes' Theorem to relate the age pdf to 
the likelihood function for the data T> = {/2, n, D*}: 

f(t\V,Z)f{V\Z) = f(V\t,l)f(t\l). (33) 

We assume that the prior probability for t, f(t\X), is constant, consistent with a constant 
neutron-star birth rate over the recent past. We evaluate the likelihood of the data in the 
same manner as Eq. |23| to find 

f t (t) = f(t\Va)=K l f dDD 4 J V P) / dVPf(z M P) )f S ^-^ P) -^ DGR) ) (34) 

where J\ft is the normalization constant. 

Application to individual pulsars is straight forward. Figure [Tl] shows the pdfs f t (t) 
and corresponding cdfs for three pulsars. Also shown in the figure, as vertical lines, are the 
conventional spindown age estimates, ts- The shapes of the individual pdfs are directly 
related to the Gaussian shapes of the assumed underlying pdfs for z , Vq P ^ and By 
inspection of Eq. [34], if measurement errors are infinitesimal and the values for birth z and 
D are fixed, then f t (t) has weight where f(z , Vq ) maximizes, in particular, near the 
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origin. If acceleration is not significant then the minimum |Vg | implies = and the 
minimum z implies to ~ sinb/nw Thus, the maximum (or mode) of f t is at an age t . The 
figures show that the resultant pdf is not symmetric, falling off more slowly toward large 
t than small t because of the effect of peculiar radial velocity and the assumed Gaussian 
form: large (giving small t) has exponentially small probability of occurrence whereas 
small VP ( giving large t) has only algebraically small probability. 

Figure |TT] shows the 68% probability ranges for the ages of the 47 pulsars in our sample 
along with the median, mode, and spindown age. As demonstrated in § [4.5| , our likelihood 
analysis shows a better fit to the pulsar data for combinations of braking index and torque 
decay time that make ages of most pulsars smaller than the conventional spindown age, 
Ts = P/2P. The conventional spin down time assumes a braking index of 3, an infinite 
torque-decay time, and birth period much smaller than the present-day period. The same 
trend is apparent in Figure [TT] which shows that the majority of objects have spindown ages 



To further demonstrate that a bias exists between chronological and spindown ages, we 
inspect the ratio of chronological age to spindown age, r = t/r s . The pdf for this ratio is 



We average f r (r) over 45 objects in our sample, excluding the Crab and Vela pulsars because 
their kinematic ages are poorly constrained. Figure |12| shows the average f r , which peaks 
at r ~ 0.4 and has a fast rise at smaller r and a slow decay for larger r. The mode of this 
pdf indicates that kinematic ages are systematically smaller than conventional spindown 
ages. This conclusion is in accord with our previous discussion, where we invoked specific 
mechanisms to maximize the likelihood in our velocity analysis and which also implied 
smaller kinematic ages. The results given in this section demonstrate this conclusion more 
directly and also show that it is a general trend in the population and not an artifact 
produced by a few extreme cases. Also, the pdf for age is calculated independently of any 
particular physical process, such as braking mechanisms other than that due to magnetic 
dipole radiation and torque decay. 



In Appendix [A] we discuss constraints and uncertainties for particular objects, especially 
those with uncertain distances and large velocities. Here we comment on the possible birth 
sites and velocities of a few objects. 

The Geminga pulsar (B0630+17) was discovered as a gamma-ray source and eluded 



rs > t - 




(35) 



7. 



CONSTRAINTS ON INDIVIDUAL PULSARS 
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detection as a periodic source until 1992 (Halpern & Holt 1992). Recent parallax and 
earlier proper-motion observations (Caraveo 1993; Caraveo et al. 1996) yield fairly tight 
determinations of its distance and proper motion. Caraveo et al. discuss possible birth 
sites for this pulsar, including work by Gehrels & Chen (1993), Frisch (1993), Smith et 
al. (1994), and Cunha & Smith (1996). The latter three authors suggest birth in the event 
that has produced a supernova-like ring around AOri, a possible companion to the pulsar 
progenitor, now requiring a radial velocity of about —700 km s -1 . Our constraints on 
the radial velocity (Table ||), V r ~ ± 178 km s _1 suggests this association is unlikely, 
essentially the same argument by Carveo et al. , who state that the velocity would make an 
unlikely small angle with the line of sight (~ 11°). 

8. A VELOCITY-MAGNETIC MOMENT CORRELATION? 

Various kick mechanisms suggest the magnetic field of a young pulsar is correlated with 
its acceleration.^ Horowitz and Piekarewicz (1996; see also Bisnovatyi-Kogan 1996) have 
discussed the role of parity violation for neutrinos interacting with electrons and nucleons in 
a strong magnetic field. Several interactions (neutrino scattering from polarized nucleons, 
from polarized electrons, and in the presence of a polarized background of electrons) can 
couple the momentum flux of the neutrinos to the B field direction with a magnitude 
proportional to B, potentially leading to the recoil of the entire neutron star. A field of 
~ 10 16 G is needed to achieve the observed velocities. A more speculative suggestion of 
Kusenko & Segre (1996a) is that neutrino oscillations are effected by the strength of the 
magnetic field. If the field is sufficiently strong (~ 10 14 G according to the original proposal, 
~ 10 16 G according to revised estimates by Qian 1997), then interconversion of v e and v T 
between the two different neutrinospheres allows some of the trapped electron neutrinos to 
convert to the tau form and escape. Neutron star recoil results from anisotropic emission 
of neutrinos and the resultant momentum impulse. The resonant interconversion depends 
on the function k ■ B and the net momentum impulse could vanish for certain perfectly 
symmetric field configurations (e.g. a strong, exactly toroidal field but not a dipole field). 
Kusenko & Segre (1996b) discussed an expected correlation in their model, v oc B for long 



3 According to the "rocket theory" of Harrison & Tademaru (1975), pulsars with off-axis magnetic dipolcs 
are accelerated as they spin down. For rapidly spinning neutron stars, the final velocity is independent of 
the magnetic moment since larger magnetic moments cause a larger force that acts for a shorter time. The 
net velocity created by the Harrison- Tademaru rocket depends on the birth period and offset of the magnetic 
axis and one would not expect to find a correlation between a pulsar's velocity today and its initial magnetic 
field. However, a relation between the offset and magnetic field strength would introduce a correlation. 
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period pulsars. Birkel and Toldra (1997) considered the correlation between transverse 
velocity and the estimated magnetic field projected along the pulsar spin axis for short 
period objects. The former find evidence in favor of the hypothesis, the latter do not. 

Some previous investigations have identified significant correlations between 
perpendicular velocities and magnetic moment (e.g. Anderson & Lyne 1983; Cordes 1986, 
1987; Bailes 1989) to varying degree and others have not (Lorimer, Lyne & Anderson 1995). 
Using the derived velocities for our sample, we tested for a correlation between the 3D 
velocities of Table |5] and various combinations of P and P by cross correlating log V and 
log PP 13 . Figure O shows the correlation coefficient plotted against (3. The correlation 



maximizes at (5 = with a correlation coefficient of 0.26, which can be achieved by chance 
at the 7.8% level. This indicates that the 3D velocity is not correlated with the magnetic 
moment (/3 = 1) or with the spindown age (J3 = —1) or any other combination that involves 
a significant contribution from P. Moreover, velocities are at most weakly correlated with 
spin period. 

Since there has been some disagreement over the significance of the correlation in the 
past, we provide the following discussion. The previously analyzed samples included objects 
which have been spun up, i.e. millisecond pulsars, and which have significantly lower space 
velocities and smaller magnetic fields than the young, pulsars considered in this paper. As 
shown in paper I, (Cordes & Chernoff 1997) the nominal rms 3D velocity for MSPs ~ 80 km 
s _1 is a factor ~ 5 smaller than that of young pulsars. It seems clear that the correlations 
between V and PP previously identified were effectively measuring this difference between 
the young and MSP populations. As argued in paper I, both the low magnetic fields and 
space velocities of MSPs are related to the required binary configuration needed in order 
for accretion driven spinup to occur. Accretion quenches or rearranges the magnetic field 
of the pulsar so as to diminish its apparent dipolar magnetic field. Binaries that are not 
disrupted have low center of mass velocities. This economical interpretation suggests that 
the apparent V-PP correlation has nothing to do with the physics of core collapse that 
produces neutron stars. Bailes (1989) has suggested one scenario where a pair of NSs 
are created from a binary. The binary evolves in such a way to produce one strong-field, 
high-velocity NS and one low-field, low-velocity NS. Our interpretation and our statistical 
results on the kinematics of pulsars are consistent with this scenario but do not require the 
particular assumptions of his model. 

Although we find no evidence of correlation between v and B, we do not regard that as 
particularly strong evidence for or against the role of neutrino-mediated parity violations 
or neutrino oscillations. The magnetic field required for these schemes is ~ 10 14 — 10 16 G, 
much larger than the inferred dipolar field today. Thus, the original field must have decayed 
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by 10 2 — 10 4 and it is unclear what relationship today's residual field bears to the field that 
acted during the epoch of acceleration. 

9. PULSAR-SUPERNOVA REMNANT ASSOCIATIONS 

Frail, Goss & Whiteoak (1994) deduce a mean transverse velocity ~ 550 km s -1 for 
12 pulsars near supernova remnants, using the apparent angular offset between pulsar and 
remnant centroid, a distance estimate, and an age estimate, taken to be the conventional 
spindown age. One association in their sample (B1757 — 24 & G5.4 — 1.2) yields a transverse 
speed ~ 1800 km s -1 . A few others were excluded on the basis of their 'extreme' velocities, 
ranging from 1600 to 3600 km s -1 , including an association between PSR B1610 — 50 and 
Kes 32 proposed by Caraveo (1993). 

Our constraints on kinematic pulsar ages (§^[) are based on objects that are much 
older than the characteristic 20 kyr age of those in the Frail et al. sample. If there is no 
field decay, the most likely braking index, n rs 4.5, implies that the Frail et al. velocities 
are biased too low by about a factor of two. The mean transverse speed in the Frail et 
al. sample would then be ~ 1100 km s _1 , which is about 3 times higher than the mean 
transverse speed in our best-fit model 2.1. If there is field decay on a timescale of 3 Myr and 
a braking index n w 2.5 (as in our best model), then the Frail et al. sample will have ages 
that are slightly larger than the conventional spindown times by a factor 3/2.5 = 1.2, thus 
reducing the mean, transverse velocity in Frail et al. by about 20% to ~ 460 km s _1 , which 
is still larger than the transverse mean (~ 331 km s^ 1 ) of our velocity distribution by 26%. 

In either case, the velocities calculated from pulsar-supernova remnant associations 
are larger than the mean pulsar distribution. This is intriguing and may indicate that 
supernova remnants are produced preferentially by the explosions that yield fast kicks. 
It is clear that a statistical analysis of the significance of the pulsar-supernova remnant 
associations and a more complete analysis of the selection effects in the pulsar surveys are 
needed. For example, some of the higher velocity objects in the Frail et al. sample may be 
spurious associations. To gain some sense of the incompleteness of the pulsar survey at high 
velocity, we can assume that the 12 associations used by Frail et al. are real and estimate a 
correction to the pulsar birth rate for the selection against high-velocity objects in pulsar 
surveys as follows. Frail et al. estimate a median velocity ~ 460 km s™ 1 that we correct to 
380 km s _1 (using a braking index of 2.5). From our velocity pdf based on 47 pulsars with 
average spindown age 2.95 Myr, we find P{V± > 380 kms -1 } = 0.22. Assuming identical 
low velocity distributions, our sample is depleted of high velocity pulsars by a factor 
0.5/0.22 ~ 2.3. For all velocities the sample has been depleted by a factor 1/0.72 ~ 1.4. 
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Thus, if the Frail et al. estimates are correct, then about 40% of pulsars are missed due to 
selection against high-velocities. 

10. ASTROPHYSICAL CONSEQUENCES 

In this section we discuss a few of the astrophysical implications of the results we have 
reached. 

10.1. Evidence for Asymmetric Supernova Kicks 

The large peculiar velocities we have estimated for pulsars — and their distribution- 
are the net result of kicks from supernova explosions, including rocket-like accelerations 
(Harrison & Tademura 1975), and disrupted orbital motion in those cases where a neutron 
star, or its progenitor, have become unbound from a binary or multiple star system. 

A considerable body of evidence supports the conclusion that supernova explosions are, 
in general, asymmetric so as to impart significant velocities to neutron stars. This evidence 
includes (1) statistical arguments based on population synthesis studies which indicate that 
symmetric explosions yield too many binary radio pulsars and too small a mean speed 
(Dewey & Cordes 1987; hereafter DC87); (2) the existence of objects, like the Guitar Nebula 
pulsar (B2224+65), which have velocities larger than that of any plausible progenitor 
binary; (3) the occurrence of geodetic precession in two NS-NS binaries, which requires 
a misalignment of the spin and orbital angular momenta through explosion asymmetry 
(Cordes & Wasserman 1984; Weisberg et al. 1989; Cordes, Wasserman & Blaskiewicz 1990; 
Arzoumanian et al. 1996) suggestive of a kick ~ 100 — 300 km s _1 ; (4) orbital precession in 
the pulsar J0045-7319, probably the result of spin-orbit coupling (Lai et al. 1995; Kaspi et 
al. 1996). Assuming aligned angular momenta in the pre-supernova state, Kaspi et al. infer 
a probable minimum kick of 100 km s _1 to achieve the current misalignment of spin and 
orbital momenta. 

Quantitatively, except for the Guitar nebula pulsar, these facts each call for kicks 
>100 km s _1 . For example, DC87 concluded that ~ 90 km s _1 kicks were needed to 
account for the then-observed ~ 100 km s _1 pulsar velocity distribution and the low binary 
frequency of pulsars. Current population synthesis models of neutron star containing 
binaries predict birth rates of wide and close radio pulsar binaries and of high mass X ray 
binaries in rough agreement with observations when kicks are drawn from a Gaussian with 
dispersion ~ 450 km s _1 (Portegies Zwart & Verbunt 1996). Without the kick, far too 
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many Be-pulsar binaries, Be- and high-mass X-ray binaries and binary radio pulsars are 
found. It is important to note that the minimum kick needed to remove the wide binaries 
is less than that needed to explain the proper motions studied in this paper. Also, the only 
direct measurement of very fast ~ 1000 km s _1 objects is that of the Guitar nebula pulsar. 
Indirect evidence comes from statistical studies of proper motion like the current one and 
that of Lyne & Lorimer (1994). Less clear evidence comes from associations of pulsars with 
some supernova remnants of known age (see §||). 



10.2. High Velocities and Escape from the Galaxy 

The current results extend previous results in a number of ways. We begin with the 
high velocity end of the pdf. The pulsar sample itself includes 5 examples (11 %) whose 
inferred 3D velocities exceed 1000 km s _1 . 

Our model fitting allows us to estimate the fraction of objects that have sufficient total 
velocity to escape the Galaxy. The escape velocity is, of course, not known with precision 
and it is a strong function of location in the Galaxy. Leonard & Tremaine (1990) place a 
lower bound on the local V e of 430 km s~ 1 . Integrations of the Paczynski potential also 
indicate that 500 km s _1 is a representative estimate for V e , so we take this as a fiducial 
value. 

In Figure |] we show the pdfs for the peculiar speed and the cumulative distribution 
functions (cdfs) of speed for three models {n g .n h = 1.1,2.1 & 3.1). The poorly-fitting, 
single-component model (1.1) suggests that 50% of the known, young pulsars will escape 
the Galaxy. The best model (2.1) and slightly- worse model (3.1) yield a more modest 
~ 20%. The cdf for pulsar speeds indicates 5-25 % above 1000 km s -1 , taking into account 
uncertainties in our fitting, as indicated in Figure || 

In interpreting these results, it is necessary to reiterate that the pulsar sample is known 
to be strongly biased by selection effects, many of which arise from the fact that pulsar 
surveys are signal-to-noise limited. Many factors contribute to the net S/N for a given 
pulsar, including period dependence of the underlying pulsar beam shapes and luminosity. 
More importantly, however, is the velocity selection that occurs (Helfand & Tademaru 
1977; Cordes 1986; Lorimer & Lyne 1994): a greater number of fast pulsars will move out 
of the detection volume centered on the Earth than will move into it. Thus, pulsar surveys 
are biased against high velocity pulsars. The magnitude of this bias depends on details of 
specific pulsar surveys. Our analysis does not incorporate the detection-related part of the 
observer's selection function (see introduction) and it is not possible to correct post facto 
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the derived pdf. We are currently working on incorporating the detection-related terms 
to the selection function (Chernoff & Cordes 1997). The velocity pdf derived by Lyne & 
Lorimer (1994) has an overall shape that is similar to our double Gaussian model, but 
predicts a larger fraction of pulsars that will escape the Galaxy, > 50%. 

10.3. Number of Separate Velocity Components 

The two component distribution, Model 2.1, we have found to be the best fit in §4.4] , 
may or may not correspond to two physical subcomponents in the Galactic population of 
radio pulsars. Our original assumption was that the sum-of-Gaussians form would suffice 
to model any actual pdf, with no implication that the individual Gaussians corresponded 
to distinct physical subpopulations. However, the fact that the three component model 
has an odds ratio inferior to the two component model and, in addition, that there is a 
significant disparity between the individual velocity parameters (crvi ~ 4oy 2 ) is suggestive 
that two distinct subpopulations may be involved. At the same time, we note that we have 
explored only Gaussian models and that testing other parametric forms will be necessary 
before drawing this conclusion. 

To test for the presence of an additional low- velocity population of pulsars, we explored 
solutions where two of the components were varied near our best fit for the two-component 
model (2.1) and a third additional component was imposed with dispersions in the range 
10 km s -1 < <Jv3 < 210 km s _1 . By varying the weights, Wi,w 2 , dispersions ayi, &V2 while 
holding the birth scale height (h zl ) fixed, we find w 3 < 0.02, 0.05, 0.1 for ov 3 < 30, 50 
and 70 km s _1 , respectively. We therefore find no evidence in the data for a low-velocity 
component that might provide a signature of symmetric supernovae in isolated pulsars or 
in wide binaries. There is only a few percent allowed for pulsars with speeds less than 50 
km s _1 , consistent with the appearance of only two objects in our sample (B1747— 28 and 
B1929+10) whose 3D speeds could be, but need not be, this small. 

If there were no intrinsic kick then isolated supernovae would produce low velocity 
neutron stars with velocities of order their progenitors' speeds, ~ 30 km s _1 . The fraction 
of NS that originate from isolated (nonbinary) progenitors is not well known. In the models 
of Portegies Zwart & Verbunt (1996) composed of 40% single stars, 60% binaries, with 
LL type kicks at birth (models SS and AK) approximately 34% of all supernova occur in 
isolated progenitors, 66% in binaries. Since our model implies that there are relatively few 
low velocity pulsars, a natural interpretation is that all neutron stars receive a substantial 
kick at birth. 
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This conclusion is at variance with the recent proposal of Iben and Tutukov (1996) that 
there are no intrinsic kicks during supernova collapse, that instead, the pulsar velocity is 
produced entirely by the breakup of binaries. To avoid overproducing low velocity pulsars, 
Iben and Tutukov posit that pulsars can arise only from stars in which mass transfer or 
tidal interactions function to spin up the star prior to collapse. Comparison with the 
Portegies Zwart & Verbunt results shows that this assumption is not sufficient to zero out 
the low velocity part of the pdf. Roughly half of the supernovae in binaries involve mergers 
of stars prior to the explosion. We anticipate that the merged objects will resemble single 
stars in their center of mass motion but will contain substantial spin angular momentum. 
If these objects form neutron stars as the population synthesis models predict and if the 
neutron stars become visible pulsars but do not receive birth kicks then they will contribute 
weight to the low velocity part of the pdf. This is not what we observe. Thus, whether or 
not single stars give rise to pulsars it appears that kicks ^50 km s" 1 are generally required 
for consistency with the pdf. 

A supernova kick of ~ 100 km s _1 has generally been called for in population synthesis 
models (DC87): otherwise too many neutron star containing binaries are created. It is 
satisfying that the binary abundance gives results consistent with the low velocity end of 
the pulsar pdf. These abundances are an intrinsically less clear cut constraint because of 
the many complexities of binary evolution including the uncertainties implicit in the initial 
conditions. On the other hand, the small size of our current sample limits the precision of 
the direct determination of the pdf at low velocity. Future proper motion measurements 
should improve the situation. 

Recent searches for old, accreting, isolated neutron stars have turned up few prospects 
(see Wang 1996 for a review and Manning, Jeffries & Willmore 1996 for limits on the 
ROSAT Wide Field Camera survey). Our results imply that fewer than 10% of neutron 
stars move with velocities less than 70 km s^ 1 (see above). This differs from the estimates of 
~ 25 — 30% used by Madau & Blaes (1994) which were based on the young pulsar velocity 
distribution of Narayan & Ostriker (1990), diffusively heated by gravitational interactions. 
Since Bondi-Hoyle accretion scales as v ~ 3 , the new velocity distribution should play an 
important role in determining the number of detected objects in X ray surveys and may 
help reconcile the paucity of sources and the theoretical predictions. The accretion rate 
onto the star depends also on the strength of the magnetic field and the rotation rate, so 
our estimates of the braking index and torque decay timescale are relevant too. 

For supernovae in binaries that do not involve mergers (roughly 30%) the resultant 
velocity distribution represents a convolution of the orbital velocity distribution with the 
intrinsic kick velocity distribution. The maximum velocity attainable with 90 km s _1 
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kicks for initial masses up to 2OM was estimated by DC87 to be about 600 km s _1 (by 
reaction [Nl], with such a small separation that the supernova shock completely destroys 
the companion) but this constituted < 0.2% of all formed neutron stars. A more significant 
8% contribution of ~ 500 km s _1 was found in reaction [N3] in which non-conservative 
evolution creates an unbound neutron star and white dwarf. These numbers are, of course, 
only examples based on the input quantities used by DC87; especially important is their 
attribution of final (pre-explosion) orbital periods after the spiral-in processes that influence 
the Nl and N3 paths; these periods were only guessed at, implying that the typical 500-600 
km s -1 velocities estimated for these paths are not hard numbers. 

Iben and Tutukov (1996) argued that the highest orbital velocity of neutron star are 
created when the primary-descended neutron star in a binary is released as a consequence of 
a second supernova explosion. They envision a semi-detached binary consisting of a 16M 
He core and an neutron star with a 1.7i? Q semi-major axis. The orbital velocity is ~ 1300 
km s' 1 and the first neutron star is ejected at ~ 900 km s _1 . Given that the initial mass 
of the main sequence progenitor of the He core is about 40M Q , close to the maximum mass 
for a star to yield a neutron star, their example provides an upper bound to the possible 
orbital velocity. Our result that 5 of the pulsars have velocities in excess of 1000 km s _1 
then provides decisive evidence that kicks of some sort are required. 

In the future, once the observational selection effects are more fully accounted for, 
it will be important to analyze the velocity pdf in detail. At the moment, we make a 
few general comments. An outstanding issue is whether the two Gaussian description 
corresponds to a single distribution with an extended high velocity tail or whether it more 
closely corresponds to two separate populations of pulsars. Our two component results 
imply low and high velocity pulsars with comparable numerical weights (based on the crude 
estimate of the selection against high velocity pulsars above). If the contribution from the 
kick and from the orbital disruption of binaries were randomly combined (in the sense of 
two random velocities drawn from individual Gaussians) then the resulting distribution 
would be unimodal in speed with a dispersion that is the root mean square of the individual 
dispersions. The current result is not unimodal. 

This may be an indication that the two known velocity sources, orbital motion and 
intrinsic supernova kick size, are correlated. For example, binaries in which super giants 
explode will generally be larger than those in which He cores explode (since mass transfer 
has occurred in the second but not the first case). If the kick size of the former is smaller 
than the latter, then the pdf for velocity may have two distinct components. In fact, 
supernovae will occur in a variety of modified stellar types in binaries (He core, CO core, 
Wolf-Rayet and super giant according to Portegies Zwart & Verbunt 1996) and more 
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complicated correlations with orbital motion can be expected. Another possibility is that 
the merged stars with no orbital motion but large angular momentum may experience kicks 
of a different characteristic size than regular or stripped stars. While all these possibilities 
are speculative, the accurate determination of the pdf, whether or not it is unimodal, 
promises to be of value in elucidating the complexities of binary evolution. 



10.4. Implications for core collapse 

An important aspect of this work is that we establish that a significant fraction of 
the pulsar velocity distribution lies above 1000 km s _1 . This should provide important 
constraints on numerical simulations of supernova explosions. 

In a supernova collapse the energy scales are 10 53 ergs for neutrino emission and neutron 
star binding energy, 10 51 ergs for mass motions plus optical display, and 1.4 x 10 49 ergs for 
a 1.4 solar mass neutron star moving at 1000 km s -1 . Thus, the energy involved in NS 
motion is but a small fraction of the total. Moreover, the NS velocity is quite small relative 
to free fall velocities in the vicinity of the core. It is interesting to note that the distance 
traveled by the core moving at a typical speed of 1000 km s _1 for an interval corresponding 
to the delay before explosion (0.07-0.42 s for SN1987a as estimated by Bethe 1993) is 
substantial compared to the characteristic dimensions of the system. Bethe estimates 
three important lengthscales: the neutrinosphere where the optical depth to neutrinos 
becomes of order unity (20 km), the gain radius where heating by the radiated neutrinos 
exceeds the gas' cooling (~ 140 km) and the free nucleon radius where alpha particles are 
completely dissociated into neutrons and protons (~ 210 km). As the protoneutron star 
forms, the accretion shock stands off at ~ 300 km and begins to accelerate outward once 
the ram pressure of the accreting material diminishes sufficiently The neutrino heating 
and convection behind the shock are essential for providing the pressure that launches the 
shock. The core's motion is slow and the material between the neutron star and the shock 
is subsonic, so whatever acceleration mechanism functions, it should make little difference 
for estimating the state of the normal matter within the shock. 

We begin by describing the limits to the generation of core motion. If the kick Av ns 

is ultimately due to the momentum impulse associated with anisotropic neutrino emission, 

then a firm limit may be expressed in terms of the total binding energy of the neutron star, 

e-g- ^ 

Av ns < e v — — -c (36) 
M ns c z 

where e v is the anisotropy in the neutrino momentum and AE is the total energy emitted, 
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which must be less than the binding energy of the final neutron star 0.15M ns c 2 . To give 
kicks of 1500 km s" 1 , we require that e u > 0.03. From numerical simulations, Janka & 
Miiller estimate that the energy that is available is significantly less than the binding energy 
- more closely the accretion of the mantle during the convective period: AE « 0.015M ns c 2 . 
Then a substantially larger asymmetry e v > 0.3 is required to give 1500 km s _1 . This 
degree of asymmetry is inconsistent with simulations. 

If the kick is due to the momentum recoil from either accreted or expelled matter, the 
momentum transfer is of order the escape velocity near the gain radius times the mass 
contained between the gain region and the accretion shock front times an asymmetry factor, 
e.g. 

AM . . 

Av ns < e ej ——v ej (37) 

where e e j is the anisotropy in the momentum carried off, AM is the mass and v e j is the 
escape velocity. From numerical simulations, Janka & Miiller estimate from the most 
extreme cases that AM/M ns < 0.2, v e j < 0.1c and e e j < 0.08 and infer Av ns < 500 km s _1 . 
Once again, a substantially larger asymmetry e e j > 0.24 would be required to give 1500 km 
s- 1 . 

These general arguments are not without loopholes, of course. They depend on the 
validity of the existing numerical simulations to set reasonable limits for e v or e e j. In all 
calculations to date the neutron star is fixed at the center of symmetry. As mentioned 
above, because the fastest neutron star motions are subsonic these approximations should 
be reasonable for gauging the size of purely hydrodynamic quantities like AM, AE, v e j and 
e e j that appear above. The fact that neutrino transport is treated with limited accuracy 
and that current numerical modeling is not done in full 3D are much more difficult issues 
to evaluate. 

As many others have discussed, it would appear that the likely physical sources for a 
kick are (1) the large neutrino flux beyond the neutrinosphere and (2) the outer envelope 
of infalling material which is not in sonic communication with the matter surrounding the 
core. Goldreich, Lai & Sahrling (1997) have suggested that core g modes in presupernova 
stars may be overstable. As the matter collapses, irregularities in the infalling envelope may 
grow to large amplitude and ultimately impart a kick by gravitational torques. Another 
possibility is envisioned by Burrows and Hayes (1996) in which the inhomogeneities of the 
infalling matter channel the supernova blast into a jet-like outflow that causes the core to 
recoil. 

It is important to point out, also, that the lack of correlation between the velocity and 
magnetic field in today's pulsars (§§) does not necessarily rule out a neutrino-magnetic field 
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interaction ultimately responsible for the recoil. If field decay has occurred, as suggested by 
our best fitting models with torque decay, then the correlation between magnetic field and 
velocity could be weakened. 



11. SEARCHING FOR THE HIGHEST- VELOCITY PULSARS 

It is almost certain that the velocity pdf we have derived underestimates the fraction 
of pulsars with V <: 10 3 km s _1 . Searching for high velocity pulsars at high galactic 
latitudes is a clear, though nontrivial, way to improve the situation. The difficulty arises 
because (1) such pulsars will be faint owing to their presumably large distances; (2) proper 
motion measurements will require VLBI techniques; given pulsar faintness, the VLBA will 
have to be augmented by the large-aperture Arecibo Telescope and Green Bank Telescope; 
(3) high-velocity pulsars will be above the thickest free electron layer in the Galaxy and 
therefore their distances cannot be estimated from the TC (or any other) model for the 
electron density. Difficulties (1) and (2) can be surmounted through suitable allotment of 
telescope time. Estimating pulsar distances is much more difficult to tackle, for inverting 
dispersion measures will never be possible for pulsars with \z\ ^ 1 kpc. Parallax observations 
with VLBI or pulse-timing techniques are also unlikely. The only possible recourse is 
development of an alternative distance scale. We are investigating the possibility that radio 
pulsars have beam components whose luminosities are determined by P and P and no other 
parameters. This 'standard candle' approach also requires that beam shapes and angular 
diameters also be predictable, in which case orientation angles between the pulsar spin and 
magnetic axes and the observer's direction can be determined, along with the predicted flux 
and, hence, distance. Work on this beam modeling will be reported elsewhere. Should beam 
modeling fail due to (e.g.) pulsar 'weather' effects, such as severe distortions of magnetic 
fields from dipolar forms by accretion or other history dependent activity, then analysis of 
high latitude pulsars must rely on statistics of large numbers of objects rather than velocity 
estimates of individual objects. 



12. SUMMARY AND CONCLUSIONS 

In this paper we have used the data on 47 pulsars to infer their velocity and z-altitude 
distributions at birth. Our likelihood method allows estimation of radial as well as 
transverse velocities, contingent on estimates for the ages and birth z's of individual objects. 
We emphasize that we have not corrected the results for the strong selection effects in the 
original pulsar surveys that have provided the sample of objects. The strongest effect on 
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our conclusions is that we have underestimated the fraction of pulsars with large velocities, 
e.g. those with V > 1000 km s _1 . We are currently developing methods, similar to those 
in this paper and in Paper I, that take selection effects into account. 

With these provisos, we have found that: 



1. All pulsars in the sample are consistent with birth within 0.3 kpc of the galactic plane; 
any apparent motion of pulsars toward the galactic plane is easily explainable as a 
combination of projection effects associated with pulsars that are relatively nearby, a 
conclusion made long ago, on a smaller sample, by Helfand & Tademaru (1977). 

2. The best fit model we have found consists of two Gaussian components in velocity, 
with dispersions ~ 175 and 700 km s _1 , combined with a birth scale height of 0.13 
kpc; we have not investigated mathematical forms other than the multicomponent 
Gausssians. 

3. With our fit, we predict that ~ 20% of pulsars like those in our sample will escape 
the Galaxy; a rough estimate indicates that this fraction may be underestimated by a 
factor ~ 2 due to selection effects. 

4. A braking index n ~ 2.5 is favored if the pulsar torque constant decays on a time 
scale ~ 3 Myr. If no torque decay occurs, then most braking indices are in the range 
3.5-6, though such a model yields a slightly inferior likelihood than does the model 
with torque decay. 

5. The constraints on braking index and torque decay time arise from the role that the 
kinematic age plays in our estimates of radial velocities. The most general statement 
is that pulsar ages are typically a factor ~ 2 smaller than the conventional spindown 
ages calculated as P/2P, which assumes birth at a period much smaller than the 
present day period and a braking index of 3, with no torque decay. 

6. We find no significant correlation between velocity and magnetic field estimates. 
Reconciling this result with previous studies is simple: previous work included strong 
field and weak field objects that had been spun up by accretion. We found in Paper 
I that millisecond pulsars are generally a low velocity population. The correlation 
found in previous work appears to signify only the different evolutionary paths that 
have led to strong field objects and spun-up, low field objects. 

7. Pulsar velocities require some sort of rocket mechanism that augments velocities that 
would obtain under the condition of symmetric supernovae in binary and nonbinary 
progenitors. No strong discrimination exists in our analysis for proposed rocket effects 
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that take place during or shortly after supernova explosions. However, the estimated 
velocities require a more efficient acceleration mechanism than has been calculated 
thus far. 

8. We suggest that the best means for pushing the limits and testing models is to search 
for, and measure the proper motions of, pulsars at high galactic latitudes. 
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APPENDIX 



A. Individual Objects 

Distances, if misestimated by the TC model, are more likely to be overestimated than 
underestimated because HII regions or other enhancements of the free-electron density will 
perturb the dispersion measure (DM) upwards. In some cases, the distance estimate can be 
grossly overestimated. 



B0531+21: The Crab pulsar clearly was born near its present z distance from the 
galactic plane (~ —0.2 kpc). However, its progenitor star may have originated much 
nearer the plane. Gott et al. (1970) argued that the Crab's progenitor and the 3.75 s 
pulsar B0525+21 were unbound from the explosion that produced the latter pulsar; 
the point of origin was the Gem I OB association. 

B0540+23: A modest V r ^ < 500 km s~ x results if the pulsar was born at z ~ —0.3 
kpc. 

B0630+17: The Geminga pulsar is one of the few with a parallax as well as proper 
motion measurement. Its speed is modest and the likely range of radial velocities 
renders unlikely any association with the star AOri, which would require a radial 
velocity of about —700 km s _1 (Caraveo et al. 1996). A histogram of Monte Carlo 
values of radial velocity (cf. §|]) shows only 0.3% of counts have radial velocity 
magnitudes as large as 700 km s _1 . 

B0736— 40: The TC model yields only a lower bound on the distance (11.0 kpc) 
because its DM is too large for its latitude. Using this lower bound, the required radial 
velocity is |V^ p )| > 10 4 km s" 1 . The distance is likely to be much smaller because an 
HII region may account for much of DM. Frail (1990) estimates the distance to be 
only 0.4 kpc based on an analysis of HII regions along the line of sight. 

B1508+55: The catalog uncertainty in distance is 60%. The transverse speed 
exceeds Vj> p ^ by a factor of 7. 

B1642— 03: Prentice & ter Haar (1969) first proposed that the DM to this pulsar is 
contributed to significantly by an HII region. Its distance is then 0.1-0.2 kpc rather 
than the 1.6 kpc from the TC model. We have used this lower distance range. 

B1706— 16: Frail (1989) argues that an HII region can account for some of the DM 
for this pulsar, which might place the pulsar significantly nearer (~ 0.1 kpc) than its 
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nominal distance, 1.27 kpc. This possibility is less certain than that for B1642— 03, so 
we have used the larger distance. 

B1933+16: No auxiliary data contradict the adopted distance or the large transverse 
speed. The radial velocity need not be large if the birth altitude ~ 0.15 kpc. 

B1953+50: The radial velocity exceeds the transverse speed by a factor ~ 3.5. No 
auxiliary data contradict the assumed distance. To reduce |V^ p )| to less than 10 3 
or 500 km s" 1 , the distance would have to be reduced < 1.1 kpc or < 0.54 kpc. 
Alternatively, for a distance equal to the midpoint, 0.5(D L + Djj) = 1.84 kpc, the age 
must be reduced by a factor ~ 4 — 6 to reduce the radial velocity magnitude to 1000 
or 500 km s _1 . 

B2154+40: The radial velocity magnitude is about twice the transverse speed. No 
auxiliary data exist to contradict the derived distance and velocities. 

B2021+51: The parallax from VLBI measurements (Campbell et al. 1996) yields a 
distance consistent with that from the TC model. The spindown age agrees with the 
kinematic age to better than a factor of two. 

B2224+65: The Guitar-Nebula pulsar is the largest, convincing pulsar velocity in 
excess of 800 km s _1 (Cordes, Romani & Lundgren 1993). At the TC distance, the 
implied radial velocity and birth z are both modest. Evidently, most of the pulsar's 
speed is transverse to the line of sight and parallel to the galactic plane. It will escape 
the Galaxy. 
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Table 1: Pulsar Proper Motion Data 



Name 


t 


b 






logTs 


D L 


Du 


Kei 


(J2000) 


(B1950) 


(deg) 


(deg) 


( ill CIS \ 

\ yr ' 


/ Tn cis \ 
\ yr / 


(yr) 


(kpc) 


(kpc) 




0139 + 5814 


0136 + 57 


129.2 


-4.0 


-11.0 + 5.0 


-19.0 + 5.0 


5.6 


2.22 


3.76 


1 


0332 + 5434 


0329 + 54 


145.0 


-1.2 


17.0 + 1.0 


-13.0 + 1.0 


6.7 


1.10 


1.86 


1 


0358 + 5413 


0355 + 54 


148.2 


0.8 


15.0 + 10.0 


-6.0 + 10.0 


5.8 


1.59 


2.69 


2 


0454 + 5543 


0450 + 55 


152.6 


7.5 


52.0 + 6.0 


-17.0 + 2.0 


6.4 


0.61 


1.03 


1 


0502 + 4654 


0458 + 46 


160.4 


3.1 


-8.0 + 3.0 


8.0 + 5.0 


6.3 


1.37 


2.31 


1 


0528 + 2200 


0525 + 21 


-176.1 


-6.9 


-20.0 ± 19.0 


7.0 ± 9.0 


6.2 


1.75 


2.95 


1 


0534 + 2200 


0531 + 21 


-175.4 


-5.8 


-16.0 + 11.0 


-2.0 + 8.0 


3.1 


1.54 


2.60 


3 


0543 + 2329 


0540 + 23 


-175.6 


-3.3 


19.0 + 7.0 


12.0 + 8.0 


5.4 


2.72 


4.60 


1 


0614 + 2229 


0611 + 22 


-171.2 


2.4 


-4.0 + 5.0 


-3.0 + 7.0 


5.0 


3.63 


6.14 


1 


0629 + 2415 


0626 + 24 


-171.2 


6.2 


-7.0 + 12.0 


2.0 + 12.0 


6.6 


3.59 


6.07 


1 


0630 - 2834 


0628 - 28 


-123.0 


-16.8 


-5.0 + 18.0 


-17.0 + 26.0 


6.4 


1.65 


2.78 


2 


0633 + 1746 


0630 + 17 


-164.9 


4.3 


138.0 + 4.0 


97.0 + 4.0 


5.5 


0.12 


0.22 


4 


0653 + 8051 


0643 + 80 


133.2 


26.8 


19.0 + 3.0 


-1.0 + 3.0 


6.7 


2.32 


3.93 


1 


0659 + 1414 


0656 + 14 


-158.9 


8.3 


64.0 + 11.0 


-28.0 + 7.0 


5.0 


0.58 


0.99 


5 


0738 - 4042 


0736 - 40 


-105.8 


-9.2 


-80.0 + 13.0 


-3.0 + 43.0 


6.6 


0.40 


0.80 


2 


0742 - 2822 


0740 - 28 


-116.2 


-2.4 


-29.0 + 0.9 


-0.1+0.4 


5.2 


1.40 


7.70 


6 


0820 - 1350 


0818-13 


-124.1 


12.6 


9.0 + 10.0 


-31.0 + 7.0 


7.0 


1.89 


3.20 


1 


0826 + 2637 


0823 + 26 


-163.0 


31.7 


61.0 + 3.0 


-90.0 + 2.0 


6.7 


0.29 


0.49 


3 


0835 - 4510 


0833 - 45 


-96.4 


-2.8 


-48.0 + 2.0 


34.9 ± 1.0 


4.0 


0.39 


0.65 


6 


0837 + 0610 


0834 + 06 


-140.3 


26.3 


2.0 + 5.0 


51.0 + 3.0 


6.5 


0.55 


0.94 


3 
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Table 1: Pulsar Proper Motion Data (continued) 



Name 


I 


1 

b 




1*6 


1 

logT 5 


D L 


Du 


Ref 


(J2000) 


(B1950) 


(deg) 


(deg) 


( vnas \ 
\ yr ' 


( vticls \ 
\ yr ' 


(yr) 


(kpc) 


(kpc) 




0908 - 1739 


0906 - 17 


-113.9 


19.9 


27.0 + 11.0 


-40.0 ± 11.0 


7.0 


0.48 


0.82 


1 


0946 + 0951 


0943 + 10 


-134.6 


43.1 


-38.0 + 19.0 


-21.0 ± 12.0 


6.7 


0.75 


1.27 


3 


1136 + 1551 


1133 + 16 


-118.1 


69.2 


-102.0 + 5.0 


357.0 + 3.0 


6.7 


0.21 


0.35 


3 


1453 - 6413 


1449 - 64 


-44.3 


-4.4 


-16.4+1.1 


-21.3 + 0.8 


6.0 


2.00 


5.00 


6 


1509 + 5531 


1508 + 55 


91.3 


52.3 


-73.0 + 4.0 


-68.0 + 3.0 


6.4 


1.49 


2.51 


3 


1559 - 4438 


1556 - 44 


-25.5 


6.4 


11.0 + 17.0 


20.0 + 50.0 


6.6 


1.50 


2.50 


2 


1645 - 0317 


1642 - 03 


14.1 


26.1 


41.0 + 17.0 


-25.0 ± 11.0 


6.5 


0.10 


0.50 


3 


1709 - 1640 


1706 - 16 


5.8 


13.7 


75.0 ± 20.0 


147.0 ± 50.0 


6.2 


0.98 


1.65 


2 


1752 - 2806 


1749 - 28 


1.5 


-1.0 


-5.0 + 17.0 


8.0 + 15.0 


6.0 


1.18 


1.99 


2 


1820 - 0427 


1818-04 


25.5 


4.7 


3.0 + 3.0 


27.0 + 3.0 


6.2 


1.61 


2.73 


3 


1825 - 0935 


1822 - 09 


21.4 


1.3 


10.0 + 19.0 


-23.0 ± 19.0 


5.4 


0.78 


1.31 


2 


1844 + 1454 


1842 + 14 


45.6 


8.1 


-9.0 + 10.0 


45.0 + 6.0 


6.5 


1.70 


2.87 


1 


1913 - 0440 


1911-04 


31.3 


-7.1 


7.0 + 13.0 


-5.0 + 9.0 


6.5 


2.48 


4.19 


1 


1919 + 0021 


1917 + 00 


36.5 


-6.2 


-2.0 + 30.0 


-1.0 ± 10.0 


6.4 


2.56 


4.33 


1 


1932 + 1059 


1929 + 10 


47.4 


-3.9 


79.0 + 6.0 


39.0 + 4.0 


6.5 


0.13 


0.22 


3 


1935 + 1616 


1933 + 16 


52.4 


-2.1 


2.0 + 3.0 


-25.0 + 5.0 


6.0 


6.11 


10.32 


7 


1948 + 3540 


1946 + 35 


70.7 


5.0 


-9.0 + 7.0 


-4.0 + 8.0 


6.2 


6.05 


10.22 


1 


1955 + 5059 


1953 + 50 


84.8 


11.6 


-23.0 + 5.0 


54.0 + 5.0 


6.8 


1.37 


2.31 


1 


2022 + 2854 


2020 + 28 


68.9 


-4.7 


-9.0 + 3.0 


-13.0 + 2.0 


6.5 


1.00 


1.69 


3 


2022 + 5154 


2021 + 51 


87.9 


8.4 


-8.1 + 0.2 


13.4 + 0.2 


6.4 


0.94 


1.59 


8 



Table 1: 



Pulsar Proper Motion Data (continued) 



Name 


i 


b 




A*<s 


logr 5 


D L 


D v 


Ref 


(J2000) 


(B1950) 


(deg) 


(deg) 


/ mas \ 
\ yr ' 


/ mas \ 

\ y r ) 


(yr) 


(kpc) 


(kpc) 




2048 - 1616 


2045 - 16 


30.5 


-33.1 


85.0 + 16.0 


-43.0 + 17.0 


6.5 


0.49 


0.83 


2 


2055 + 3630 


2053 + 36 


79.1 


-5.6 


-3.0 + 7.0 


3.0 + 3.0 


7.0 


4.28 


7.23 


1 


2113 + 2754 


2110 + 27 


75.0 


-14.0 


-23.0 + 2.0 


-54.0 + 3.0 


6.9 


1.05 


1.78 


1 


2157 + 4017 


2154 + 40 


90.5 


-11.3 


18.0 + 1.0 


-3.0 + 1.0 


6.8 


4.25 


7.19 


1 


2219 + 4754 


2217 + 47 


98.4 


-7.6 


-12.0 + 8.0 


-30.0 + 6.0 


6.5 


1.89 


3.19 


7 


2225 + 6535 


2224 + 65 


108.6 


6.8 


144.0 + 3.0 


112.0 + 3.0 


6.0 


1.54 


2.60 


1 


2354 + 6155 


2351 + 61 


116.2 


-0.2 


22.0 + 3.0 


6.0 + 2.0 


6.0 


2.53 


4.28 


1 



References: 1. Harrison et al. (1992) 2. Fomalont et al. (1992) 3. Lyne, Anderson & Salter 
(1982) 4. Caraveo et al. (1996) 5. Thompson & Cordova (1994) 6. Bailes et al. (1990) 7. 
Downs & Reichley (1983 8. Campbell et al. (1996) 
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Table 2: Parameter Search Ranges 



parameter units min max range 

wi 1 Awi = 1 

w 2 1 Aw 2 = 1 

h zl kpc 0.01 0.51 Ah zl = 0.5 

h z2 kpc 0.01 0.51 Ah z2 = 0.5 

a Vl km s- 1 25 2025 Acr Vl = 2000 

o V2 km s- 1 200 2000 Aa V2 = 1800 

a V3 km s- 1 100 400 Aa V3 = 300 
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Table 3: 3D Velocity PDF Models 



Model w 1 


w 2 




h Z 2 


0Vi 


<7V2 


°V3 


■Nparms 


log£L 


Odds 


n g .n h 




(kpc) 


(kpc) 


(km s- 1 ) 


(km s- 1 ) 


(km s _1 ) 








t s < 10 Myr: 


Npsr ~ 


= 47 
















1.1 1 




0.13 




300 






2 


-263.50 


1 


2.1 0.8 


0.2 


0.13 




175 


700 




4 


-255.92 


1Q 6.3 


2.2 0.82 


0.18 


0.13 


0.10 


175 


750 




5 


-255.89 


1Q 5.5 


3.1 0.3 


0.15 


0.13 




156 


750 


200 


6 


-255.90 


1Q 5.4 


t s < 5 Myr: N psr = 


38 
















1.1 1 




0.13 




300 






2 


-208.10 


1 


2.1 0.87 


0.13 


0.13 




175 


700 




4 


-201.94 


10 51 


t s < 1 Myr: N psr = 


12 
















1.1 1 




0.13 




200 






2 


-55.99 


1 


2.1 0.8 


0.2 


0.13 




175 


400 




4 


-55.95 


0.36 
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Table 4: Best-fit Parameters for Model 2.1 

parameter value units 

0.8318:8 - 
0.131831 k P c 

175±lg km s" 1 
700l?£ km s" 1 



h zl 
(Jvi 

&V2 



Table 5: 



Pulsar Velocity Components 
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Table 5: 



Pulsar Velocity Components (continued) 
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Table 5: Pulsar Velocity Components (continued) 
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Fig. 1. — Z- velocity vs. time for objects moving near the Sun in the Paczynski potential. 
Curves are shown only for V z > 0. 

Fig. 2. — Lines showing most likely values for and zq. For three pulsars (B0540+23, 
B 1449-64 and B2224+65), dashed lines indicate the 68% confidence interval. 

Fig. 3. — Most likely values for V} p ^ and zq for four pulsars using three different values of 
braking index. Three lines are shown for each value of braking index to designate the peak 
likelihood and the 68% confidence regions. Heavy solid lines: n — 4. Light solid lines: n — 3. 
Long dashed lines: n = 2.5. Dotted lines: n = 2. 

Fig. 4. — Contour plot of the log likelihood function for the single Gaussian model 
(n g .nh = 1.1) as a function of the rms velocity and scale height. Contours are spaced 
by log 2. 

Fig. 5. — Contour plots of the log likelihood function for the double Gaussian model 
{rig.rih = 2.1) as a function of values for pairs of parameters and for 2D slices through 
the best-fit model. Contours are spaced by log 2. 

Fig. 6. — Marginalized probability density functions for the four parameters of the double 
Gaussian model, n g .rih = 2.1. 

Fig. 7.— One-dimensional cuts through the likelihood functions for a double-Gaussian 
velocity pdf where the torque decay time tk is allowed to vary for different values of braking 
index. The heavy solid line for n = 2.5 yields the maximum likelihood solution at r K ~ 3 
Myr. The two heavy dashed lines, for n = 4 and 4.5, yield the maximum likelihood if there 
is no decay, tk — > oo. 

Fig. 8. — (Top) Cumulative distribution functions for the velocity magnitude at birth using 
best fit parameter values for models 1.1 (dotted line), 2.1 (heavy solid line), and 3.1 (dashed 
line). The light solid lines represent model 2.1 evaluated with parameter values that are 
±1<7 from the best fit values. The vertical dashed line marks 500 km s _1 , the nominal speed 
of escape from the Galaxy at the solar circle. (Bottom) Differential probability density 
functions for models 1.1, 2.1 and 3.1 using best fit parameter values . 
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Fig. 9. — Scatter plots of radial velocity against birth altitude, zo, and perpendicular 
speed, V±, for three pulsars. The values result from the Monte Carlo solutions derived using 
Eq. |32|. (a) B1953+50; (b) B2154+40; (c) B2224+65 



Fig. 10. — Probability density functions and cumulative distribution functions for the ages 
of three pulsars. Vertical lines indicate the conventional spindown ages obtained from Eq. |2] 
with a braking index, n = 3. 



Fig. 11. — Plot showing the acceptable range of possible ages for each of 47 pulsars, calculated 
using the pdf of Eq. |34|. The horizontal line indicates the interval enclosing 68% of the 
probability. The x indicates the median of each pdf, the open circle the mode of the pdf, 
and the solid circle indicates the spindown age Ts = P/2P. For the Crab and Vela pulsars, 
the kinematic age is unconstraining because of their young ages; for them, we show only 
their spindown times. 



Fig. 12. — Probability density for the ratio of pulsar age to spindown age, r = t/r s , based 
on the individual pdfs for 45 pulsars. The Crab and Vela pulsars have been excluded from 
our sample because their kinematic ages are poorly constrained. 



Fig. 13. — Plot of correlation coefficient between log V and log P a P 13 vs. (3. 
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